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We present a new type of a thermodynamically stable magnetic state at interfaces and surfaces of chiral magnets. 
The state is a soliton solution of micromagnetic equations localized in all three dimensions near a boundary and 
contains a singularity, but nevertheless has a finite energy. Both features combine to a quasi-particle state for 
which we expect unusual transport and dynamical properties. It exhibits high thermal stability and thereby can be 
considered as promising object for fundamental research and practical applications in spintronic devices. We provide 
arguments that such a state can be found in different B20-type alloys e.g. Mni_ 2 ,Fea;Ge, Mni_a;Fea:Si, Fei_ 3 ,Co 3 ,Si. 


Many branches of modern physics focus on the sys¬ 
tems, which exhibit a rich variety of metastable states. 
Special attention is drawn to systems characterized by 
high energy barrier, which are able to persist in a par¬ 
ticular metastable state for an unbounded or noticeably 
long time. Such systems are within the scope of in¬ 
terest for optical, magnetic and electronic data storage 
technology. Basically, almost all modern electronics and 
spintronics are based on handling and manipulating such 
metastable states. For instance, the magnetic domains 
in magnetic data storage devices are the best examples 
of such metastable states separated by high-energy barri¬ 
ers. In hard disk drives high energy barriers are provided 
by a high magneto-crystalline anisotropy of the mate¬ 
rial, which makes at the same time the magnetic domains 
moveless and affixed to their positions. The disadvantage 
of such data storages is that the detector has to be moved 
to the place where the data bit is stored or the whole 
storage media has to be moved towards the detector. 

Modern spintronics offers many alternative ideas for 
the next-generation magnetic data storage devices that 
are based on utilizing movable metastable states. The 
most promising candidates are domain walls in nanowires 
[1] and localised magnetic vortices [2] also known as chi¬ 
ral magnetic skyrmions [3, 4]. The latter appear in mag¬ 
nets with Dzyaloshinskii-Moriya interaction (DMI) [5, 6] 
as an isolated metastable state or as a ground state, when 
skyrmions condense into a lattice. 

From the point of view of the field theory, skyrmions 
are two-dimensional solitons known as an exception to the 
Hobart-Derrick theorem [7] due to the presence of Lifshitz 
invariants [8]. Due to their unusual transport properties, 
and compact size, which at some conditions can be re¬ 
duced up to a few nano-meters [9] or even a few atomic 
distances [10], such types of chiral votices are considered 
as promising particle for information technology as well 
as an interesting object for fundamental research in non¬ 
linear physics. 

In this paper, we present a new type of metastable 
particle-like state - a hybrid particle composed of a 
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smooth magnetization vector field and a magnetic singu¬ 
larity. It is a three-dimensionally localized soliton of the 
nonlinear equations for a unit vector field. Such a quasi¬ 
particle is more compact and it turns out to be energet¬ 
ically more favorable than the chiral magnetic skyrmion 
in a wide range of parameters. 

Our approach is based on a classical spin model of the 
simple cubic lattice described by the following Hamilto¬ 
nian [11] consisting of three contributions: the Heisenberg 
exchange interaction, the Dzyaloshinskii-Moriya interac¬ 
tion and Zeeman term 


E = -j'^ni Uj • [n^xiij] (1) 

(ij) (ij) i 


where {ij) denote summation over all nearest-neighbour 
pairs, rii = Mi//rs is a unit vector of the magnetic mo¬ 
ment at lattice site i, J is the exchange coupling constant 
and Dy is the Dzyaloshinskii-Moriya vector defined as 
D.y = Dvij, D is the DMI scalar constant and is a 
unit vector pointing from site i to site j, H is an external 
magnetic field applied along 2 -axis, perpendicular to the 
film plane. Note, we assume isotropic DMI which means 
that the value jD^- j is fixed for all three spatial directions. 
Periodical boundary condition is applied in xy-plane. 

In order to preserve the generality of the results of our 
discrete model and to be able to compare them to the ear¬ 
lier developed theory based on continuum approximation 
[2, 12, 13, 14, 15] we use the following notations 


A J 

Lo = 47r— = 2na—, Hu 


2M,A 




( 2 ) 


where Lu is the lowest period of incommensurate spin 
spiral which exists below the critical field Hu; A and T> 
are micromagnetic constants of exchange and DMI respec¬ 
tively; Ms is the magnetization of the material; a is the 
lattice constant. 

We use a nonlinear conjugate gradient method 
with adaptive stereographic projections implemented on 
NVIDIA CUDA architecture for direct minimization of 
the Hamiltonian (1)[16]. 
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Figure 1. Phase diagram of the ground state for isotropic 
helimagnetic thin film of the thickness L at magnetic field H 
applied normally. The dashed area within the conical ground 
state corresponds to the isolated skyrmion as the lowest energy 
metastable state. Arrows indicate the values of the fixed pa¬ 
rameters used in the energy calculation presented in Figs. 3a-d. 
The star symbol corresponds to the parameters used for the 
calculation of non-zero anisotropy case, see Figs. 3e, f; the 
circles - for the calculation of the energy barriers, see Fig. 4. 

In order to identify the range of parameters corre¬ 
sponding to the stable configurations we have calculated 
the magnetic phase diagram for the ground state, see 
Fig. 1. Here, solid lines correspond to the first order 
phase transition between helical state [17], skyrmion lat¬ 
tice (SkL), conical phase, and saturated state. Horizontal 
dashed line between conical and saturated states corre¬ 
sponds to the second order phase transition. Under the 
condition Ld is one order of magnitude or more larger 
than a, discrete (1) and continuum [16] approaches are 
consistent. 

Results presented in Fig. 1 are consistent with the con¬ 
tinuum theory [18] and experimental data on the direct 
observation of SkL in thin films of cubic helimagnetic al¬ 
loys as Fei_ 2 ;Coa;Si [11], Mni_a;Fea;Ge [9], Mni_a;Fea;Si 
[19] and pure MnSi [20] and FeGe [21]. The mechanism 
of SkL stabilization in thin films of cubic helimagnet dif¬ 
fers from two-dimensional systems with interface induced 
DMI [22] and anisotropic bulk helimagnets [23, 24]. Ac¬ 
cording to Ref. [18], the free boundary of the layer which 
provides an additional degree of freedom to the magnetic 
system and changes the common energy balance are re¬ 
sponsible for the SkL stabilization. In particular, due to 
the lack of exchange coupled neighbors on the surface the 
spin structure of an isolated skyrmion tube (SkT) exhibits 
a small twist of magnetization with respect to the surface 
normal, as schematically shown in Fig. 2a. Such an in¬ 
duced twist propagates through the whole thickness of the 
sample, accumulates the energy gain of DMI contribution 
along the z-axis and results in significant reduction of the 
total energy of the SkT. As a result, within a certain range 
of applied field the energy of the SkL becomes lower than 
the energy of the conical phase. The equilibrium value 
of the angle for such a surface twist is defined by the 
competition between positive contributions of exchange 
coupling and two negative contributions: i) in-plane DMI 
- between the spins in the same plane and ii) out-of-plane 
DMI - between the spins along z-axis. For the SkT and 
ChB this angle (po < 15° (see Fig. 2). Recently the effect 


of the chiral surface twist has been proven experimentally 
[25, 26]. 



Figure 2. Schematic representation of the vector field and 
cross sections of corresponding isosurfaces (riz = 0) for 
skyrmion tube (SkT) (a) and chiral bobber (ChB) (b). The 
small (green) sphere indicates the position of the chiral Bloch 
point. The spins in red marked box on the right illustrate the 
magnetization distribution in the conical phase, which sur¬ 
rounds these localized metastable states. Circular blue arrows 
indicate the sense of magnetization rotation within the conical 
phase and surface induced twist. For the exact spin structure 
of SkT and ChB see Supplementary Movies 1 and 2 [16] . 

Everywhere within the conical phase the lowest en¬ 
ergy metastable state corresponds to the that shown in 
Fig. 2b. The exception is the small dashed area in Fig. 1 
where SkT dominates. 

The vector field corresponding to such a solution is 
characterised by the shape of the isosurface = 0, which 
is close to a shape of paraboloid, see Fig. 2. An essential 
element of the corresponding magnetic structure is a sin¬ 
gularity situated at a finite distance P from the surface, 
see the green sphere at the bottom of the red isosurface. 
This type of singularity is known in ferromagnets as Bloch 
point (BP) or hedgehog [27], and in other condensed mat¬ 
ter systems e.g. ultra cold gas [28] and spin ice [29] as 
magnetic monopole. In helimagnets, due to the presence 
of DMI, the vector field around the singularity exhibits 
a chiral structure. Thus we use the name chiral BP to 
refer to it. Earlier it has been shown that some dynam¬ 
ical processes in chiral magnets are accompanied by the 
appearance of this kind of singularities [30, 31]. We suc¬ 
ceeded to get an analytical expression for the vector field 
in the close vicinity to the chiral BP, which gives an op¬ 
portunity to estimate its energy contribution to the total 
energy of the system [16]. 

In the cross section with the plane parallel to the 
layer surface its spin structure mimics the magnetization 
distribution of skyrmion with the diameter diminishing 
with the distance from the surface. We found that the 
value of the penetration depth always satisfy the condi¬ 
tion P < Ld/ 2 with an accuracy up to ±a/2. Because of 
the essential chirality of such a state and its localization 
close to the surface with the finite penetration depth, like 
a fishing bobber at a water surface, we use term chiral 
bobber (ChB) to refer to this object. 
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Figure 3. Energy dependencies for isolated skyrmion tube (SkT, open circles) and chiral bobber (ChB, red solid circles) on 
applied magnetic field for fixed film thickness (a)-(b), on film thickness for fixed applied field (c)-(d), and as a function of varying 
values of cubic (e) and uniaxial anisotropy (f) for fixed thickness and field. The energies are presented relative to the pure conical 
phase. Note, the cases in (a)-(d) correspond to the isotropic case, = 0. The top and bottom scales in (a-b) and (c-d) 

are identical but given in absolute and reduced units, respectively. For these calculations we used J = 1, D = 0.15 (Ld = 41.89a) 
and a size of the simulated domain of 256 x 256 x spins, where L„ = 40 4- 280. 


We calculated the energy dependencies of the ChB and 
SkT at applied field and varying thickness, see Figs. 3a-d. 
Here, we use the same color notation as in the phase dia¬ 
gram, Fig. 1: white and grey areas correspond to the sta¬ 
bility range of conical phase and SkL, respectively, dashed 
area denotes the range within the conical ground state 
where the SkT is lower in energy than the ChB. As fol¬ 
lows from the Figs. 3a-d almost everywhere within the 
conical phase the ChB has energy much lower than the 
SkT. The energy difference between the ChB and SkT 
becomes most pronounced by varying the thicknesses, see 
Figs. 3c, d. The energy of the SkT increases linearly with 
the thickness. At the same time, due to the localisation 
of the ChB near the surface its energy is almost thickness 
independent. It shows significant variation only at low 
fields and very small thickness, outside the range of the 
conical phase stability, see insets in Fig. 3d. 

In order to estimate the stability of the ChB in 
anisotropic systems we add to the model Hamiltonian 
Eq.(l) the following term 

Ek = -Y^{ K^nl^ + ATc } > (3) 

i 

where n^x is cc-component of unity vector n at the lattice 
site i, Kc is the cubic anisotropy constant and AT^ is the 
constant of uniaxial anisotropy induced in the system, e.g. 
by the interfaces or under external stress. In Figs. 3e, f 
the energies of the ChB and SkT are presented as func¬ 


tions of varying anisotropy constants at fixed applied field 
H = O.SHy) and thickness L = 3Ld. We use reduced units 
of anisotropy where ATq = a^I?^/(4A) = D'^/{2J). 

As follows from the figures, in the whole range of ex¬ 
istence within the conical phase the ChB remains ener¬ 
getically more favorable than the SkT. The star symbols 
in Figs. 3e and f corresponds to the collapse of the ChB. 
In Fig. 3e the narrow gray region denotes the range of 
SkL stability, which is consistent with previous theoreti¬ 
cal study [24]. 

Due to the presence of free boundaries and singulari¬ 
ties as an essential element of the solution the topology of 
such object should be described in frame of relative homo- 
topy groups [32]. Since, our boundary conditions impose 
no restrictions on the magnetization, as an order param¬ 
eter, the corresponding group is trivial, 7r2(S^,§^) = 0. 
It means that the singularity can be pushed out on the 
surface and then smoothed. It is important to empha¬ 
size that when certain restrictions to the order parameter 
are imposed at the boundary it may lead to formation of 
topologically nontrivial states, as in case of a boojum in 
A-phase in ^He [33, 34]. Despite the fact that solution 
for the ChB is non-topological, it shows the high energy 
barrier. 

In order to estimate the energy barriers we have calcu¬ 
lated the minimum energy path (MEP) with the geodesic 
nudged elastic band (GNEB) method [35], which is an ex¬ 
tension of the original nudged elastic band method [36, 37] 
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to magnetic systems. In the GNEB method, the curvature 
of the configurations space of a magnetic system arising 
from a constraint on the length of magnetic moments is 
taken into account. The path between the metastable 
states is presented by the discrete sequence of system 
snapshots called images. These images are then brought 
to the MEP in the curved configuration space. 



Figure 4. Minimum energy path calculated for Li = 2.07 Ld, 
and L 2 = 2.79Z/D at H/Hy> = 0.8, Kc = Ku = 0. Reaction 
coordinate is an order parameter which has a meaning of the 
relative distance between the states in the configuration space 
[35]. The local minima a, c, and e correspond to the energies 
of SkT, pair of ChBs, and single ChB, respectively, see also 
schematic representation in the corresponding insets. The dif¬ 
ference between the energy values in the local minimum and 
nearest maximum corresponding to the saddle points (see 6, 
d, and /) defines the energy barrier responsible for the sta¬ 
bility of corresponding state. The reference level {E = 0) is 
the energy of the global minimum g corresponding to the pure 
conical phase [16]. 

Examples of MEPs are presented in Fig. 4. For this 
calculations we used the following parameters: J = 1, 
D — 0.45 (Ld = 13.96a) and domain size of 30 x 30 x 30 
and 30 X 30 X 40. The elimination of the SkT occurs via 
nucleation of a pair of chiral BPs (point h) moving to¬ 
wards opposite surfaces and results in appearance of two 
repulsive ChBs (point c) which have energy higher than 
two isolated ChB. 

As follows from Fig. 4, the energy barriers for the ChB 
and the SkT have comparable height, AEef = 0.43AEab 
for Li and AEef = 0.28AEab for L 2 . From that we 
may conclude that the thermal stability range for ChB 
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Appendices 

Al. Micromagnetic analysis for point singularities in isotropic chiral mag¬ 
nets 

In this section we derive an analytical expression for the vector field in the close vicinity to the point singularity which 
allowed us to identify the role of Dzyaloshinskii-Moriya interaction in its energetics. 

It is easy to show that the model Hamiltonian Eq. (1) of the main text can be reduced in continuum approximation 
to the following functional [40, 41]: 

£ = [ W£{r)dr, W£[r) = A{dxT0? + dyVL^ + dzTo?)+V [n-voin) + HMs{l - Uz) (Al.l) 

Jv 

where n = n(r) is a continuous unit vector field defined everywhere except at the singular points. A and D are 
micromagnetic constants for exchange and Dzyaloshinskii-Moriya interactions, respectively, H is the applied magnetic 
field, and is the magnetization of the material (the total magnetic dipole moment per unit volume). For the case of 
a simple cubic lattice micromagnetic constants are related to the microscopic constants of the discrete model as 

J/(2a),D = D/a^Ms =/is/a^ (Al.2) 

where a is the lattice constant, fis is the magnetic dipole moment of the magnetic atom. 

In spherical coordinate parametrisation, the components of the unit vector field are defined as 

n = (sin(0) cos($), sin(0) sin($), cos(0)), 

where 0 = 0(r) and $ = $(r). 


5 



To minimize the functional Eq. (Al.l) one needs to solve corresponding Euler-Lagrange equations which in this 
case has a following form: 


'a (sin(20)(V$)2 - 2A0) - 2Psin(0) (n • V$) + HMs sin(0) = 0, 
-A (2 cos(0)(V0 ■ V$) + sin(0)A$) +V(n- V0) = 0. 


(A1.3) 


Let’s assume that the singular point coincide with the origin of spherical coordinate system (r, i9, <p). Then, in the 
close vicinity to this point, the Taylor expansion of 0(r) and $(r) reads 


°° /T) \ 

/V \ 


(A1.4) 

(A1.5) 


In the assumption that magnetization distribution around the singular point has an axially symmetric structure, the 
first three terms in Eqs.(A1.4, A1.5), which satisfy the Eqs. (Al.3) are 


and 


to = 2 arctan (ccot(^?/2)), 
sin((/3i) 


U — —- 


■ sin(i9), 


^ (1-c2)cos(vJi)2 . ,,, (2 c+(1 + c2)cos((/3i)) , 

t 2 = - — -sm(«) — - ---^ sm(2i>)+ 


16c 

/ c(l - c^) cAHMs 

Vl2(l + c2) + 


32c(l + A) 
cot(i9/2) 

6T>2 J 1 + A cot(t?/2)2 ’ 


(A1.6) 


/o 

h 

/2 


(1 — c^) cos(:/9i) 2c + (1 + c^) cos((^i) 


cos('d), 


4c 4c 

sin (9 (c® + c^ + c^ + l) cos {ipi) + 6c® — 8c® + 6c) 
96c2 (c2 + 1) 

(l - c^) sin {ipi) ((c^ + l) cos [ifi) + c) 


8c2 


cos(i5)+ 


(-1 - c^) sin {(fii) ((c2 + l) cos {ipi) + 2c) 

32c2 (c2 + 1) 


cos(2'd). 


(A1.7) 


where c and (pi are arbitrary free constants, which are not defined by the asymptotic expansion but they can be 
determined numerically. Note, the assumption of an axisymmetric structure, strictly speaking, is satisfied only at 
magnetic field above the saturation of the conical phase, H > Hu = /{2MsA). 

Using the Eqs. (A1.4 - A1.7) we have got an expression for the energy contribution from the volume inside the 
sphere of radius R surrounding the singularity: 


R TT 27r ^ 

£s{R) = J J J r^sm{'&)dipd'ddr = 8nAR + /3j7® + 0{R'^), 

0 0 0 


(A1.8) 


where 

(1 + c^ + 2c • cos(v?i))(l — c^ + 4 c 2 ln(c)) 4A77Msc2(l — c^ + 2 ln(c)) 

P= (l-c2)3 P2(l_c2)2 ■ 

In particular, for H — Hu, the domain of definition for jd is —1.06 < /3 < 1.32, but the exact value can be derived only 
by numerical methods. Here it is important to note, that the constants c and pi, in fact, are strictly dependent on the 
ratio H/Hu¬ 
la the limiting case H = 0 and 7? —>■ 0, the expression in Eq.(A1.8) shows a linear dependence on R and converges 
to the energy dependence derived earlier for singularity in pure ferromagnets [42]. Thus, the energy contribution to the 
total energy related to the small volume around the singularity resulting from the Dzyaloshinskii-Moriya interaction 
is proportional to i?®. Erom above one may conclude that the energetics of singular points in chiral magnets is more 
intricate than in pure ferromagnets. 
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A2. Direct energy minimization technique 

For the most general case discussed in the main text the discrete model Hamiltonian is 


E= Eq- 

(d> (u> 


[ui X rij] rii -K^ 




-K, 


E 


■ nly + nl^) 


(A2.1) 


with the assumed constraint 


2 I 2 


-ni,z = 


(A2.2) 


which means that the order parameter space is a two-dimensional sphere and each vector is considered as a point on 
it, see Fig. A2.1. For the case of such constraint a few different minimization technique have been proposed [43, 44, 45]. 
Here we use nonlinear conjugate gradient (NCG) method with Polak-Ribiere formula for updating of the conjugate 
direction. To conserve the constrain (A2.2) it is convenient to use, instead of n^, ny,nz, the stereographic projection 
defined by the pair of corresponding coordinates 71 , 72 , see Fig. A2.1 . Such a stereographic projection is defined for 
each point on the sphere, except at the sphere’s poles from which the points are projected onto the target plane. The 
projection of the point with 0 —?> 0 (0 —?> tt) is not defined for the projection from the north pole (from the south pole) 
because \/ 7 i -f 7 I —t 00 . This is quite inconvenient for a computer implementation where the values of the variables are 
always restricted to a finite range. Here we propose an alternative approach based on adaptive stereographic projections 
which allows to overcome the problem of the poles inertia [43]. Contrary to one pole stereographic projection, we 
consider both type of projections: from the north pole 


27 I 


(N) 


^ + iiry + {iry 

and from the south pole: 


271' 


iS) 


27; 


(N) 




i + iiYr + hYy 


Uv = 


27f^ 


l + (7W)2 + (7yv;) 


l + (7p^)^ + (7f^)^^ 




W^2 


(N). 


Uz = 


1-(7P^)^-(7F^)^ 


l + (7p^)2. 




(A2.3) 


(A2.4) 



Figure A2.1. Stereographic projection from the unit sphere onto the equatorial plane 2 = 0. Point P on the sphere can be 
projected from the north pole N to point P' or alternatively from the south pole S to point P" onto the target plane. 

Corresponding coordinate homeomorphisms (inverse formulas) read: 


\ {N} 


7i^^ = » 2 x/(l - riz) 
72^^ = %/(l - riz) 


(A2.5) 


\ {S'} ^ 


7^^ = ^•x/(l + Uz) 
=%/(l + nz) 


(A2.6) 


From the topology point of view, such a set of two coordinate charts is an atlas for . In our implementation of the 
NCG method, we allow an independent switching between the two types of projections for each spin, see Supplementary 
movie 4 [16]. In particular, for initial configuration we use north pole projection Eqs. (A2.3) for points with initial value 
riz < 0 and south pole projection Eqs. (A2.4) for those which have > 0, see Fig. A2.2a. On each NCG iteration 
step each spin accepts new coordinates with a small increment with respect to the previous state. Some spins which 
are initially situated in the southern semisphere can enter the northern semisphere and vice versa, see Fig. A2.2b. At 
some iteration one or few points may appear too far away from its initial semisphere, see Fig. A2.2c. The criterion for 
that can be e.g. \/ 7 i + 7 I > Fc, where Fc is some arbitrary chosen value. In our implementation we choose it such 
that \/ 7 i + 7 I = Fc which corresponds to = 0.6 for the spins defined by the north pole and = — 0.6 for south 
pole projection respectively. As soon as such event occurs at least for one spin, then all the points are switched to other 
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set of variables in the same way as at initial step, see Fig. A2.2d. Thus some coordinates ( 71 ^^ 72 ^^) are changed into 
corresponding ( 7 ^'^^, 72 '^^) and vice versa: 


(AT) 

7i = 


(S) 

7i = 


7i 


(S) 


(7P^)^ + (7r)^’ 

7 !-) 


(AT) 

72 = 


(S') 

72 = 


(S) 

72 


(7p¥ + (7f^)^’ 

7f^ 


(A2.7) 


(A2.8) 


hrr + H T hrr + iirr 

Then the NCG starts again with the new set of stereographic coordinates 71 , 72 . Note that during the minimization 
process such a global switching occures very seldom or never, it always depend on the initial configuration. 





Figure A2.2. Illustration for the adaptive algorithm of the switching between stereographic projections. Big blue dots are spins 
defined by stereographic projection from the north pole, small red dots are spins defined by projection from the south pole, see 
also Fig. A2.1. Dotted curves define critical distances from the pole of projection, see definition for Fc in the text, (a) the initial 
state, (b) the state after several iterations using the NCG method, (c) the state at which one marked point has exceeded the 
critical distance, (d) the same state as (c) but with redefined variables of the stereographic projections, see Eq. (A2.7, A2.8). 

In order to distinguish ( 71 ^^ 72 ^^) from ( 7 P\ 72 ^^) without introducing additional logical variable which define 
the type of corresponding projection for each spin e.g. N (true) and S (false), we used the following trick. In our 
stereographic projection scheme the center of the unit sphere is assumed to be located at some arbitrary position (xq, 
j/o) such that \/Xq + i/q > Fc and xq > 0, j/o > 0. Thereby, the sign of the variables 71 ,72 stored in the computer 
memory can be used as the indicator for the type of the projection, while in Eqs.(A2.7,A2.8) the absolute values of the 
variables are used and a shift of the basis (see a;o,i/o in Fig.A2.1) has to be taken into account. 

Such an adaptive stereographic projection scheme i) allows one to obtain supreme performance of the NCG method, 
ii) fully satisfies the constraint Eq. (A2.2), and iii) allows to avoid the problem with the inertia of the poles [43]. 
Contrary to other methods, see e.g. Refs. [44, 45], another important advantage of our approach is the absence of any 
free scalar parameter which would require an additional strategy to replace this free parameter by the optimal value to 
archive the best performance. 

The algorithm has been implemented in a NVIDIA CUBA architecture for massive parallel calculations and real time 
visualization and has been used for the calculations on the graphic cards with microprocessors GKI04 and GKllOB. 
In order to achieve maximum performance we used single-precision floating-point format and implemented Kahan 
summation and optimized parallel reduction algorithm [46] in order to reduce numerical errors. 

As initial conditions we used cone phase with vortex-like cylindrical inset: 


©ini = TT _|_ y^/Lo^ , 

<i)ini = arctan(j//a:) -b 7r/2 


{^/x^ + < Ld) 


(A2.9) 


where Ld = 27raJ/jDy j, see also Eq.(2) in the main text. For the calculation of the chiral bobber with the NGG 
method, the height of the cylinder has been restricted by 0.65Ld, while the thickness L of the simulated magnetic film 
has been varied in the range between ILd and 7 Ld. 


A3. Monte-Carlo simulations 

In order to identify the probability for spontaneous nucleation of the chiral bobber driven by ordinary thermal fluctu¬ 
ations we also perform simulated annealing [48] with a Monte Garlo method. We implemented the classical Metropolis 
algorithm [47]. In the simulation of the annealing process we use an exponential rule for the temperature sequence: 

r,+i=0.95T,. (A3.I) 

In the Supplementary Movie 3 [16] we present real-time simulations of the annealing process, performed for the following 
set of parameters: J = 1, D = 0.15, ^sH = 0.018 and domain size of 128 x 256 x 100 spins. Detailed description of the 
movie is presented below. 
















The simulation starts with a random spin distribution and an initial temperature of Tq = SJ/Ub, above the Curie 
temperature for the used set of system parameters. The temperature was gradually reduced down to T — 0.05 J/Zcb and 
then the NCG minimizer has been used instead of the Monte Carlo iterations in order to complete the full relaxation 
process faster. It is done only for the reason of better visualization. Performing the full relaxation with Monte Carlo 
requires much more iterations than the relaxation with NCG minimizer. 

At T = 0.6J/kB on the top surface on the right hand side of the simulated domain one can already see a few bright 
spots corresponding to the nucleated chiral bobbers, which is illustrated by scanning through the thickness. After the 
full relaxation the spin structure shows the conical state with three chiral bobbers inside, two on the top and one on 
the bottom surface. Finally we show iso-surfaces with Hz = const, which illustrate the localization and the relative 
positions of the chiral bobbers. 

The number of Monte Carlo steps as well as the variation in the domain size and coupling constants obviously affects 
the probability for the chiral bobber nucleation. As far as we pursuit the purpose to find the regime for metastable 
state nucleation we used relatively fast cooling with a reduced number of Monte Carlo steps as compared to the slow 
cooling regime when the purpose is to find the ground state. With the set of parameters mentioned above we have 
found the optimal number of Monte Carlo steps to be 75 per spin at each temperature step which allows to observe 
the spontaneous nucleation of the chiral bobber with the relatively high probability. However, even for this fast cooling 
regime the spontaneous nucleation of metastable skyrmion tube has not been observed in our simulations. 


A4. Geodesic nudged elastic band method 






Figure A4.1. The minimum energy path (MEP) identical to presented in the main text. Fig. 4. Insets show images corresponding 
to the extrema on the MEP marked with a — g. (a) isolated skyrmion tube, (b) nucleation of two chiral Bloch points, (c) 
metastable pair of chiral bobbers, (d) collapse of chiral bobber located on bottom surface via pushing out the chiral Bloch point 
to the bottom surface, (e) the single chiral bobber, (f) collapse of the single chiral bobber, (g) the conical ground state. Arrows 
indicate approximate position of the chiral Bloch point. 

In the calculation of the minimum energy path (MEP) with the geodesic nudged elastic band (GNEB) method 
we follow exactly the algorithm in Ref. [49]. We started with the relaxed isolated skyrmion and the pure conical 
state which are used as initial and final configurations, respectively, for the desired discrete MEP. We used a linear 
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interpolation with small stochastic noise to set intermediate states and construct an initial path between the initial and 
final configurations. 

The typical number of intermediated states (or images) was set to 15. Instead of the velocity projection optimization 
on a sphere [49] we used another method [50] which we found to be more convenient for our implementation, but possibly 
less efficient. In particular, for each path optimization step of each spin at the lattice site i of image k, we solve the 
following partial differential equation 

= -arifcy X [rife,, x Dfc,i], (A4.1) 

which is equivalent to the Landau-Lifshitz-Gilbert equation containing only the damping term, but with the recession 
term totally neglected. The 3A^-dimensional vector Dfe = (Dfc,i, Dfc, 2 , Dfc,Ar) (TV is the number of spins) in Eq. 
(A4.1) controls the relaxation of each image towards the MEP and is defined as 

Dfe = (Heff,fc(^^fc) (Heff,fe(nfe) • ffe)Tfe) + Ffe, (A4.2) 

where Heff,fe, Ffe, and ffe are SA^-dimensional vectors of the effective field, the so-called spring force, and the normalized 
tangent vector laying in the tangent space of image k, respectively. For detail see Ref. [49]. Following the algorithm 
Ref. [49] we used path splitting between nearest metastable states and relaxed each of local minima independently. In 
Fig.A4.1 we present the same MEP as in the main text together with exact view of the images corresponding to extrema 
on the MEP. 
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